#Leamos el archivo con la ayuda de geopandas:
import os
os.environ['USE_PYGEOS'] = '0' # only if you have pygeos

import geopandas as gpd
from  fiona import listlayers

#Guardamos los mapas
MexicoMaps='https://github.com/AyeliAlp/world_mexico/raw/main/maps/worldMaps.gpkg'

#Vemos las capas que hay en el mapa
listlayers(MexicoMaps)

#Nos aseguramos de que tengan el mismo CRS
countries.crs.to_epsg()==cities.crs.to_epsg()==rivers.crs.to_epsg()

#Creación de subconjuntos de mapas para mantener el país (Mexico):
#Guardamos al pais Mexico
mexico=countries[countries.COUNTRY=='Mexico']

#Recortamos
mexico_rivers = gpd.clip(gdf=rivers,mask=mexico)
mexico_cities= gpd.clip(gdf=cities,mask=mexico)


#Graficamos 
base = mexico.plot(facecolor="lightgreen")
mexico_rivers.plot(edgecolor='blue', linewidth=0.95,ax=base)
mexico_cities.plot(marker='+', color='red', markersize=15,ax=base)

#La unidad está en grados:
mexico.crs.axis_info

#Reproyectamos el mapa 
mexico.to_crs(3587).crs.axis_info

#Las operaciones que requieren CRS proyectado funcionarán:
#Calculamos el centroide del mapa proyectado
mexico.to_crs(3587).centroid

#Trazamos la version proyectada

base3587=mexico.to_crs(3587).plot()
mexico.to_crs(3587).centroid.plot(color='red',ax=base3587)

#Usamos el crs 3587 como opción de emergencia para reproyectar un mapa. Sin embargo, 
#para una opción más precisa, es mejor buscar las preparadas explícitamente para una ubicación particular del mundo. 
#Puede solicitar un crs por país en la pagina indicada en paideia: 
#recommended for Mexico (metros)
mexico.to_crs(6366).crs.axis_info

#Al calcular el centroide de la version 6366 nos damos cuenta que este no tiene error
mexico.to_crs(6366).centroid

#Reproyectamos:

base6362=mexico.to_crs(6366).plot()
mexico.to_crs(6366).centroid.plot(color='red',ax=base6362)

#Sabemos que mexico_6366 es un multipoligono
mexico_6366

#Solo necesitamos los bordes del 
mexico_6366.boundary

#Trazamos solo el limite del mapa proyectado
mexico_6366.boundary.plot() 

#Comprobamos el tipo de datos
type(mexico_6366.boundary)

#Algunas operaciones en geopandas requieren GDF o GS. Si necesita un GDF en lugar de un GS:
#Convertimos en GDF
mexico_6366.boundary.to_frame()

#Obtenemos un GDF muy simple, y si queremos agregar algo de información:
#Convertimos
mexico_6366DF=mexico_6366.boundary.to_frame() 

#Nueva columna
mexico_6366DF['name']='Mexico' 

#Cambiamos el nombre de la columna de geometría
mexico_6366DF.rename(columns={0:'geometry'},inplace=True) 

#establecemos la geometría (el nombre no es suficiente)
brasil_5641DF = brasil_5641DF.set_geometry("geometría")
mexico_6366DF = mexico_6366DF.set_geometry("geometry")

# Verificamos:
mexico_6366DF.crs

mexico_6366DF

#Notamos esto: (Explicar profesor)
mexico_states.crs, mexico_municipalities.crs

#Sin embargo, no tienen información de crs:
#Trazamos: 
fig, (ax1, ax2) = pyplot.subplots(ncols=2, sharex=False, sharey=False, figsize=(12,12))

mexico_states.plot(ax=ax1, facecolor='lightgrey', edgecolor='black')
mexico_municipalities.plot(ax=ax2, facecolor='lightgrey', edgecolor='black',linewidth=0.2)

mexico_states

#Como estamos usando el crs 6366, lo utilizamos tambien para los estados de mexico
mexico_states.to_crs(6366)

#Conocemos la proyeccion real de la geomretria de los estados
mexico_states.geometry

#Dissolving
#Hemos creado un subconjunto del mapa usando clip (subconjunto con respecto a los límites de otro mapa). También puede crear subconjuntos usando loc o iloc (de pandas); 
#permítanme ahora subconjunto usando cx de geopandas para hacer subconjuntos por ubicación:

centroMexico=mexico_6366.centroid


fig, axs = pyplot.subplots(nrows=2,ncols=2,sharex=False, sharey=False,figsize=(8,8))

axs[0,0].set_title('Above_Centroid, LeftOf_Centroid',fontsize=5)
mexico_municipalities_6366.cx[:centroMexico.x,centroMexico.y:].plot(ax=axs[0,0]) 

axs[0,1].set_title('Above_Centroid, RightOf_Centroid',fontsize=5)
mexico_municipalities_6366.cx[centroMexico.x:,centroMexico.y:].plot(ax=axs[0,1]) 

axs[1,0].set_title('Below_Centroid, LeftOf_Centroid',fontsize=5)
mexico_municipalities_6366.cx[:centroMexico.x,:centroMexico.y].plot(ax=axs[1,0]) 

axs[1,1].set_title('Below_Centroid, RightOf_Centroid',fontsize=5)
mexico_municipalities_6366.cx[centroMexico.x:,:centroMexico.y].plot(ax=axs[1,1])

#Creamos un poligono con una de estas partes
#En lugar de:
mexico_municipalities_6366.cx[centroMexico.x:,:centroMexico.y]

#Creamos
Mexico_topLeft=mexico_municipalities_6366.cx[centroMexico.x:,:centroMexico.y].dissolve()
#Y obtenemos
Mexico_topLeft

#Aparte de la columna de geometría, podemos observar que el resto está mal, entonces:

Mexico_topLeft['zone']='NW'
Mexico_topLeft=Mexico_topLeft.loc[:,['ADM0_ES','zone','geometry']]   
Mexico_topLeft

#Este es el trazo

Mexico_topLeft.plot()

#Tenemos que tener en cuenta que con la funcion CX es posible que no obtengamos una particion perfecta:

rowsCX=(mexico_municipalities_6366.cx[:centroMexico.x,centroMexico.y:].shape[0],\
mexico_municipalities_6366.cx[centroMexico.x,centroMexico.y:].shape[0],\
mexico_municipalities_6366.cx[:centroMexico.x,centroMexico.y].shape[0],\
mexico_municipalities_6366.cx[centroMexico.x,centroMexico.y].shape[0])
sum(rowsCX)

#Entonces, tienes polígonos duplicados en las zonas (Preguntar al profesor):

sum(rowsCX)-mexico_municipalities_6366.shape[0]

#Creación de datos espaciales
#Obtendrá Líneas y Polígonos como mapas con seguridad, pero ese puede no ser el caso con los puntos. Descargamos un archivo CSV 
#con información sobre los aeropuertos de Mexico de este sitio web, lo guardaré en mi carpeta de datos:
import pandas as pd 
infoairports=pd.read_csv(os.path.join("data","mx-airports.csv"))

#Visualizamos algunas filas:

infoairports.iloc[[0,1,2,3,-4,-3,-2,-1],:] #head and tail

#Realizamos limpieza:

#Borramos la primera fila que no nos sirve
infoairports.drop(index=0,inplace=True)
infoairports.reset_index(drop=True, inplace=True)
infoairports.head()

#Mantenemos las columnas necesarias

infoairports.columns.to_list()

#Verificamos el tipo de las variables:
infoairports.info()

#Alguno de los formatos
numericCols=['latitude_deg', 'longitude_deg','elevation_ft']
infoairports[numericCols]=infoairports.loc[:,numericCols].apply(lambda x:pd.to_numeric(x))

#Ahora
infoairports.info()

#Trazamos el mapa con los aeropuertos

base = mexico.plot(color='white', edgecolor='black') #unprojected

infoairports.plot.scatter(x = 'longitude_deg', y = 'latitude_deg',ax=base)

#Ahora trazamos

base = mexico.plot(color='white', edgecolor='black')
airports.plot(ax=base)

#Recordamos el tipo:
type(airports), type(infoairports)

#Como funciona trazamos
airports.to_crs(6366).plot()

#Verificamos que tenemos diferentes tipos de aeropuerto
airports_6366['type'].value_counts() # this will not work: airports.type.value_counts()

#Graficamos de una manera dinamica
from folium import LayerControl


m = airports_6366[airports_6366.kind=='small_airport'].explore(color="red",name="small",show=False)
m = airports_6366[airports_6366.kind=='medium_airport'].explore(m=m, color="blue",name="medium",show=False)
m = airports_6366[airports_6366.kind=='large_airport'].explore(m=m, color="black",name="large",show=True)
m = airports_6366[airports_6366.kind=='seaplane_base'].explore(m=m, color="green",name="seaplane",show=False)
m = airports_6366[airports_6366.kind=='closed'].explore(m=m, color="white",name="closed",show=False)
m = airports_6366[airports_6366.kind=='heliport'].explore(m=m, color="orange",name="heliport",show=False)

LayerControl(collapsed=False).add_to(m) #optional

m

#Creando el casco convexo
#En algún momento puede que tengas la necesidad de crear un polígono que sirva de envolvente a un conjunto de puntos. Como sobre, debe tener solo ángulos convexos. Ese sobre es el casco convexo:

# en un polígono
Mexico_topLeft

Mexico_topLeft.plot()

Mexico_topLeft.convex_hull

#Graficamos el poligono convexo 
Mexico_topLeft.convex_hull.plot()

#Observe que el casco convexo funciona para cada fila de la geometría. Brazil_topLeft era una fila de geometría, por lo que el casco convexo era sencillo.

#Esto podría ser diferente si tiene puntos:
large_airport=airports_6366[airports_6366.kind=='large_airport']
#Varias filas
large_airport

#En esta situación, como desea el casco convexo de los grandes aeropuertos, debe combinar todos los puntos:
large_airport.unary_union

#Esta es una geometría multipuntos
gpd.GeoSeries(large_airport.unary_union)

#Ahora sí obtenemos el casco convexo
large_airport.unary_union.convex_hull

#Observamos el tipo:
type(large_airport.unary_union.convex_hull)

#Pero sí es una serie
gpd.GeoSeries(large_airport.unary_union.convex_hull).plot()

#Convertimos la GS a una GDF
LargeAirport_hull=gpd.GeoSeries(large_airport.unary_union.convex_hull).to_frame()
LargeAirport_hull

#Los siguientes pasos
LargeAirport_hull.rename(columns={0:'geometry'},inplace=True)
LargeAirport_hull['name']='large airports hull' # optional
LargeAirport_hull=LargeAirport_hull.set_geometry('geometry')
LargeAirport_hull.crs="EPSG:6366"
#Verificamos:
LargeAirport_hull.crs

#Usamos el GDF en el trazado
base=mexico_6366.plot(facecolor='yellow')
large_airport.plot(ax=base)
LargeAirport_hull.plot(ax=base,facecolor='green',
                       edgecolor='white',alpha=0.4,
                       hatch='X')

LargeAirport_hull.explore()

#Comprobación de la validez de la geometría
Las #geometrías se crean de manera que pueden aparecer algunos problemas, especialmente en (múltiples) polígonos. Verifiquemos si nuestros mapas recientes sobre estados y municipios son válidos:

# no válido
#Observamos que no hay ningun dato invalido
mexico_states_6366[~mexico_states_6366.is_valid]

#Es difícil ver lo que está mal. Consigamos algo de información:

# ¿lo que está mal?

from shapely.validation import explain_validity, make_valid

explain_validity(mexico_states_6366[~mexico_states_6366.is_valid].geometry)

#Alguna invalida?
MxSt6366_valid[~MxSt6366_valid.is_valid]

#¿Qué pasa con el otro geodataframe?

mexico_municipalities_6366[~mexico_municipalities_6366.is_valid]

explain_validity(mexico_municipalities_6366[~mexico_municipalities_6366.is_valid].geometry)

#Variedades?
mexico_municipalities_6366['validity']=[x.split('[')[0] for x in mexico_municipalities_6366.geometry.apply(lambda x: explain_validity(x))]
mexico_municipalities_6366['validity'].value_counts()

#Resolviendo el problema:
mexico_municipalities_6366.drop(columns=['validity'],inplace=True)

MxMun6366_valid=mexico_municipalities_6366.copy()

MxMun6366_valid['geometry'] = [make_valid(row)  if not row.is_valid else row for row in MxMun6366_valid['geometry'] ]
#any invalid?
MxMun6366_valid[~MxMun6366_valid.is_valid]

#La solución que obtuvimos puede ayudar con algunas técnicas avanzadas, pero también puede darnos algunos problemas adicionales. Tenga en cuenta que una vez que geopandas resolvió el problema, ha creado colecciones:

[x for x in MxMun6366_valid["geometry"]]